isochrone
Yuxuan
31/05/2021
Import libraries
library(tidyverse)
library(ggplot2)
library(leaflet)
library(RColorBrewer)
library(scales)
library(lattice)
library(dplyr)
library(sf)
library(tidyverse)
# For pretty knitting
library(lemon)
knit_print.data.frame <- lemon_print
knit_print.tbl <- lemon_print
knit_print.summary <- lemon_printImport and organize data sets
## Import raw Travel Time Matrix (ttm)
ttm <- read.csv('../../data/clean/ttm.csv')
df<-read_csv('../../data/score_sets/newest_long_scores.csv')## Parsed with column specification:
## cols(
## fromId = col_double(),
## type = col_character(),
## weight = col_character(),
## nearest_n = col_character(),
## score = col_double(),
## pop = col_double()
## )# convert Ids from double to factor
ttm$fromId <- as.factor(ttm$fromId)
ttm$toId <- as.factor(ttm$toId)
target_amenities <- c('gallery', 'museum', 'library or archives', 'theatre/performance and concert hall')
amenities <- read.csv('../../data/clean/vancouver_facilities_2.csv') %>% filter(type %in% target_amenities)
# clean
amenities <- amenities[,c(1,4)] # only need id and type columns
amenities$id <- as.factor(amenities$id) # convert to factor
amenities$type <- as.factor(amenities$type) # convert to factor
# Import weight
dest_wts <- read.csv('../../data/amenity_weight/poi_index.csv')
ttm <- ttm %>% inner_join(amenities, by = c('toId' = 'id'))
head(ttm)## fromId toId avg_unique_time sd_unique_time type
## 1 59150004004 10 99.76316 5.364721 museum
## 2 59150004004 15 72.48718 3.401794 museum
## 3 59150004004 157 96.69231 3.001349 gallery
## 4 59150004004 1759 106.82051 4.388213 museum
## 5 59150004004 1760 46.58974 2.642944 gallery
## 6 59150004004 1822 76.64103 3.990035 museum### only consider the nearest amenity
nearest_1_ttm <- ttm %>%
group_by(fromId, type) %>%
summarise(avg_time = min(avg_unique_time),
toId=toId[which.min(avg_unique_time)],
sd_time = sd_unique_time[which.min(avg_unique_time)])## `summarise()` has grouped output by 'fromId'. You can override using the `.groups` argument.nearest_1_ttm<-na.omit(nearest_1_ttm)
nearest_1_ttm%>%filter(type=='museum')%>%group_by(toId)%>%filter(avg_time<15)## # A tibble: 1,984 x 5
## # Groups: toId [85]
## fromId type avg_time toId sd_time
## <fct> <fct> <dbl> <fct> <dbl>
## 1 59150028019 museum 14.5 9570 1.71
## 2 59150039004 museum 12 9569 0.427
## 3 59150039005 museum 13 9569 0.427
## 4 59150039006 museum 14 9569 0.427
## 5 59150040003 museum 11 9570 0.427
## 6 59150040004 museum 6 9570 0.427
## 7 59150040005 museum 10 9570 0.427
## 8 59150040006 museum 11 9569 0.427
## 9 59150040009 museum 12 9570 0.427
## 10 59150041004 museum 5 9569 0.427
## # … with 1,974 more rowsAssign level to each group
Group1: 0-15
Group2: 15-30 Group3: 30-45 Group4: 45-60
# check the histogram of avg time
plot(hist(nearest_1_ttm$avg_time))
# create a new column useing ifelse
# between function is inclusive <=
nearest_1_ttm$time_group= with(nearest_1_ttm, ifelse(avg_time<15,"<15",
ifelse(between(avg_time,15,29) , "15-30",
ifelse(between(avg_time,30,44),"30-45","45-60"))))
nearest_1_ttm$time_group<-as.factor(nearest_1_ttm$time_group)
nearest_1_ttm## # A tibble: 57,282 x 6
## # Groups: fromId [14,343]
## fromId type avg_time toId sd_time time_group
## <fct> <fct> <dbl> <fct> <dbl> <fct>
## 1 59150004004 gallery 35.8 3569 1.90 30-45
## 2 59150004004 library or archives 35.7 9567 1.72 30-45
## 3 59150004004 museum 38.5 9570 1.41 30-45
## 4 59150004004 theatre/performance and concer… 48.2 7239 2.60 45-60
## 5 59150004005 gallery 37.9 3569 1.60 30-45
## 6 59150004005 library or archives 37.5 9567 1.52 30-45
## 7 59150004005 museum 40.4 9570 1.44 30-45
## 8 59150004005 theatre/performance and concer… 50.8 7239 2.21 45-60
## 9 59150004006 gallery 38.9 3569 1.51 30-45
## 10 59150004006 library or archives 38.5 9567 1.54 30-45
## # … with 57,272 more rows# nearest 1 export
write_csv(nearest_1_ttm,"../../data/score_sets/nearest_1_ttm.csv")# Import shape file data
# keep necessary columns and rows for shp
canada_dbs <- st_read("../../code/Visualizations/census2016_DBS_shp/DB_Van_CMA/DB_Van_CMA.shp", stringsAsFactors = FALSE)## Reading layer `DB_Van_CMA' from data source `/Users/yuxuancui/Desktop/MDS/data599/w2020-data599-capstone-projects-statistics-canada-transit/code/Visualizations/census2016_DBS_shp/DB_Van_CMA/DB_Van_CMA.shp' using driver `ESRI Shapefile'
## Simple feature collection with 15197 features and 27 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 4001643 ymin: 1957237 xmax: 4068119 ymax: 2032663
## projected CRS: PCS_Lambert_Conformal_Conicvan_dbs <- data.frame(canada_dbs[which(canada_dbs$CMANAME == "Vancouver"), ])
van_dbs$ID <- seq.int(nrow(van_dbs))
clean_van_dbs <- van_dbs[, c(1, 28, 29)]
head(clean_van_dbs)## DBUID geometry ID
## 1 59150244005 MULTIPOLYGON (((4031013 200... 1
## 2 59150244008 MULTIPOLYGON (((4031720 200... 2
## 3 59150372006 MULTIPOLYGON (((4019999 200... 3
## 4 59150372007 MULTIPOLYGON (((4019845 200... 4
## 5 59150373001 MULTIPOLYGON (((4019658 200... 5
## 6 59150373002 MULTIPOLYGON (((4019485 200... 6# join data into a single dataframe
# convert back to sf object
van_dbs_scores <- left_join(clean_van_dbs,nearest_1_ttm , by = c('DBUID' = 'fromId'))
van_dbs_scores_sf <- st_as_sf(van_dbs_scores)
van_dbs_scores_st <- st_transform(van_dbs_scores_sf,crs = 4326)
van_dbs_scores_st## Simple feature collection with 58136 features and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -123.4319 ymin: 49.00192 xmax: -122.4086 ymax: 49.57428
## geographic CRS: WGS 84
## First 10 features:
## DBUID ID type avg_time toId sd_time
## 1 59150244005 1 gallery 35.71795 8216 4.3222110
## 2 59150244005 1 library or archives 13.00000 7270 0.4268328
## 3 59150244005 1 museum 54.53846 5738 3.9725780
## 4 59150244005 1 theatre/performance and concert hall 36.07692 3080 4.3551820
## 5 59150244008 2 <NA> NA <NA> NA
## 6 59150372006 3 gallery 19.07692 9580 1.0854200
## 7 59150372006 3 library or archives 8.00000 5134 0.4268328
## 8 59150372006 3 museum 18.66667 4783 1.0596260
## 9 59150372006 3 theatre/performance and concert hall 19.07692 9581 1.0854200
## 10 59150372007 4 gallery 19.28205 9580 1.0748000
## time_group geometry
## 1 30-45 MULTIPOLYGON (((-122.9625 4...
## 2 <15 MULTIPOLYGON (((-122.9625 4...
## 3 45-60 MULTIPOLYGON (((-122.9625 4...
## 4 30-45 MULTIPOLYGON (((-122.9625 4...
## 5 <NA> MULTIPOLYGON (((-122.9683 4...
## 6 15-30 MULTIPOLYGON (((-123.078 49...
## 7 <15 MULTIPOLYGON (((-123.078 49...
## 8 15-30 MULTIPOLYGON (((-123.078 49...
## 9 15-30 MULTIPOLYGON (((-123.078 49...
## 10 15-30 MULTIPOLYGON (((-123.0799 4...polyg_subset <- van_dbs_scores_st%>%filter(type=="museum")
score_vec <- polyg_subset$time_group
# colour palette
pal_fun <- colorFactor(
palette = c("#e30606", "#fd8d3c", "#ffe669", "#cdff5e"),
levels = unique(polyg_subset$time_group)
)
# interactive popup
p_popup <- paste0("<strong>Accessibility:</strong>", score_vec)
# Add add-ons to maps
leaflet(data =polyg_subset) %>%
addPolygons(
stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun(time_group) , # set fill colour with pallette fxn from aboc
fillOpacity = 0.4, smoothFactor = 0.2, # aesthetics
popup = p_popup) %>% # add message popup to each block
addTiles() %>%
setView(lng = -122.8, lat = 49.2, zoom = 11) %>%
addLegend("bottomleft", # location
pal=pal_fun, # palette function
values=~score_vec, # value to be passed to palette function
title = "Vancouver museum Transit Accessibility") # legend title